

function [EEG]=cvc_physio(chan,EEG)

% CVC_PHYSIO
%
%  Synopsis
%  ========
%
%                 
%  [EEG] = cvc_physio(chan,EEG,ALLEEG,CURRENTSET)
%  
%
%  -- Author: Jeppe Christensen --
%     Center for Visual Cognition, University of Copenhagen.
%     2012 - June 2012
%
%  Purpose
%  =======
%   
%   Estimating mean heart rate from 10 beats from plethysmographic recording.
%   Adding channel to dataset channel n + 1 and scaling the plethysmographic recording   
%
%  Inputs
%  ======
%
%  chan - dataset channel holding plethysmograpic recording.
%  EEG - EEGLAB dataset with plethysmographic recording.
%  ALLEEG - EEGLAB dataset holding information of all EEG sets.
%  CURRENTSET - Current dataset to be updated with heart rate. 
%  
%  Outputs
%  =======
%
%   EEG, EEG struct output to workspace updating
%
%  With inspiration from
%  ==========
%	peakdetect.m
%	Tom McMurray
%	mcmurray@teamcmi.com
fprintf('OBS!! Run "eeglab redraw" to update the dataset current active in EEGLAB gui.');


signal=-1*EEG.data(chan,:);
signal=smooth(signal,5);
fs=EEG.srate;

if ~nargin<4
   if isempty(signal)
      pospeakind=[];
      negpeakind=[];
      return
   end
else
    return
end

sizsig=size(signal);

%	while signal is unsupported, enter supported signal or return for empty outputs

while isempty(signal)|~isnumeric(signal)|~all(all(isfinite(signal)))...
      |length(sizsig)>2|min(sizsig)~=1
   signal=input(['signal is empty, nonnumeric, nonfinite, or nonvector:\nenter '...
         'finite vector or return for empty outputs\n']);
   if isempty(signal)
      pospeakind=[];
      negpeakind=[];
      return
   end
   sizsig=size(signal);
end

sizsig1=sizsig(1);
lensig=sizsig1;

lensig1=lensig-1;
lensig2=lensig1-1;

%	if signal length is 2, return max/min as positive/negative polarity peaks

if ~lensig2
   [sig,pospeakind]=max(signal);
   [sig,negpeakind]=min(signal);
   disp('2 element signal graph suppressed')
   return
end

%	generate difference signal

difsig=diff(signal);

%	generate vectors corresponding to positive slope indices

dsgt0=difsig>0;
dsgt00=dsgt0(1:lensig2);
dsgt01=dsgt0(2:lensig1);

%	generate vectors corresponding to negative slope indices

dslt0=difsig<0;
dslt00=dslt0(1:lensig2);
dslt01=dslt0(2:lensig1);

%	generate vectors corresponding to constant intervals

dseq0=difsig==0;
dseq01=dseq0(2:lensig1);
clear difsig

%	positive to negative slope adjacencies define positive polarity peaks

pospeakind=find(dsgt00&dslt01)+1;

%	positive slope to constant interval adjacencies initiate positive polarity peaks

peakind=find(dsgt00&dseq01)+1;
lenpeakind=length(peakind);

%	determine positive polarity peak terminations

for k=1:lenpeakind
   peakindk=peakind(k);
   l=peakindk+1;
   
%	if end constant interval occurs, positive polarity peak exists
   
   if l==lensig
      pospeakind=[pospeakind;peakindk];
      
%	else l<lensig, determine next nonzero slope index
      
   else
      dseq0l=dseq0(l);
      while dseq0l&l<lensig1
         l=l+1;
         dseq0l=dseq0(l);
      end
      
%	if negative slope or end constant interval occurs, positive polarity peaks exist
      
      if dslt0(l)|dseq0l;
         pospeakind=[pospeakind;peakindk];
      end
   end
end


clear dseq01
clear dsgt0 peakind

%	if initial negative slope occurs, initial positive polarity peak exists

if dslt00(1)
   pospeakind=[1;pospeakind];
   
   
%	else initial constant interval occurs, determine next nonzero slope index

else
   k=2;
   dseq0k=dseq0(2);
   while dseq0k
      k=k+1;
      dseq0k=dseq0(k);
   end
   
%	if negative slope occurs, initial positive polarity peak exists
   
   if dslt0(k)
      pospeakind=[1;pospeakind];
   end
end
clear dsgt00 dslt0 dslt00 dseq0

%	if final positive slope occurs, final positive polarity peak exists

if dsgt01(lensig2)
   pospeakind=[pospeakind;lensig];
end
clear dsgt01 dslt01

%	if peak indices are not ascending, order peak indices

if any(diff(pospeakind)<0)
   pospeakind=sort(pospeakind);
end

%   only count peaks from positive phase

pospeakind(signal(pospeakind)<mean(signal)+std(signal))=[];

%   from inter-beat-intervals to heart rate

heart_rate=diff(pospeakind/fs);

%   filter out artifactual beats (ectopic ect.)

heart_rate(heart_rate<0.5 | heart_rate>1.6)=NaN;

%   mean heart rate over 10 beats 

heart_rate_mean=smooth(heart_rate,10);

% upsamling mean heart rate to EEG timescale

tmp(1:pospeakind(1))=nanmean(60./heart_rate_mean);
for i=1:length(pospeakind)-1,tmp(pospeakind(i)+1:pospeakind(i+1))=60./heart_rate_mean(i);end
tmp(pospeakind(end)+1:length(EEG.data(1,:)))=nanmean(60./heart_rate_mean);

% adding extra channel to EEG struct with mean heart rate data

EEG.data(chan,:)=EEG.data(chan,:)/1000;
EEG.data(end+1,:) = tmp*10;
EEG.nbchan = size(EEG.data,1);
if ~isempty(EEG.chanlocs)
	EEG.chanlocs(end+1).labels = 'Hr';
end;
[~,EEG]=eeg_store([],EEG);




